(-1.5, 1.5)
(-1.5, 1.5)
Ingeniería Biomédica
2025-05-05
Let \(x[n]\) be a discrete-time signal.
The Z-Transform is defined as:
\[X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}\]
Where: - \(z \in \mathbb{C}\) is a complex variable - \(z = re^{j\omega}\)
Causal Signals
ROC is outside outermost pole
Anti-Causal Signals
ROC is inside innermost pole
\[H(z) = 1.00 \cdot \frac{(z - 0.50)}{(z - 0.90)}\]
(-1.5, 1.5)
(-1.5, 1.5)
If the ROC includes the unit circle, \(|z| = 1\), then:
\[X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}\]
So the Fourier Transform is a special case of the Z-Transform.
Let:
\[x[n] = a^n u[n]\]
Where: - \(a \in \mathbb{R}\) - \(u[n]\) is the unit step function (0 for \(n<0\), 1 for \(n\geq 0\))
Apply the definition:
\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \sum_{n=0}^{\infty} (az^{-1})^n\]
This is a geometric series:
\[X(z) = \frac{1}{1 - az^{-1}} = \frac{z}{z - a}\]
A geometric series is a sum of terms where each term is multiplied by the same constant:
\[S = \sum_{n=0}^{\infty} r^n = 1 + r + r^2 + r^3 + \cdots\]
The value of \(r\) determines whether this sum converges (has a finite limit) or diverges.
Understand why this series:
\[\sum_{n=0}^{\infty} r^n\]
converges if and only if:
\[\boxed{|r| < 1}\]
Let’s consider the sum up to the \(N\)-th term:
\[S_N = \sum_{n=0}^{N} r^n = 1 + r + r^2 + \cdots + r^N\]
This has a closed-form expression:
\[S_N = \frac{1 - r^{N+1}}{1 - r} \quad \text{for } r \neq 1\]
To find the sum of the infinite series, take the limit as \(N \to \infty\):
\[S = \lim_{N \to \infty} S_N = \lim_{N \to \infty} \frac{1 - r^{N+1}}{1 - r}\]
If \(|r| < 1\), then:
\[r^{N+1} \to 0 \quad \text{as } N \to \infty\]
So the sum becomes:
\[S = \frac{1}{1 - r}\]
✅ The geometric series converges.
❌ In all cases: the series diverges
\[S = 1 + 0.5 + 0.25 + 0.125 + \cdots = \frac{1}{1 - 0.5} = 2\]
Every term adds less. The sum “flattens out.”
The Z-Transform often gives us geometric series like:
\[\sum_{n=0}^{\infty} (az^{-1})^n\]
This converges only if:
\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]
So, understanding convergence of geometric series = understanding ROC in Z-transforms.
| Condition | Behavior | Result |
|---|---|---|
| \(|r| < 1\) | Terms shrink | Series converges |
| \(|r| \geq 1\) | Terms grow or oscillate | Diverges |
\[\sum_{n=0}^{\infty} r^n = \frac{1}{1 - r} \quad \text{if } |r| < 1\]
For convergence of the geometric series:
\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]
Therefore, the ROC is:
\[\boxed{|z| > |a|}\]
\[x[n] = (0.5)^n u[n]\]
Z-Transform:
\[X(z) = \sum_{n=0}^{\infty} (0.5)^n z^{-n} = \frac{z}{z - 0.5}\]
Region of Convergence:
\[\boxed{|z| > 0.5}\]
| Signal | Z-Transform | ROC |
|---|---|---|
| \(x[n] = a^n u[n]\) | \(\frac{z}{z - a}\) | \(|z| > |a|\) |
| Example: \(a = 0.5\) | \(\frac{z}{z - 0.5}\) | \(|z| > 0.5\) |
Let \(x[n] = a^n u[n]\), where \(|a| < 1\)
\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \frac{1}{1 - az^{-1}}, \quad \text{ROC: } |z| > |a|\]
A difference equation relates input and output values at different time steps.
\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]
Common in: - Digital filters (FIR, IIR) - Signal models in ECG, EEG analysis - Implementation in real-time biosignal systems
The Z-Transform turns time shifts into powers of \(z^{-1}\):
| Time Domain | Z-Domain |
|---|---|
| \(x[n]\) | \(X(z)\) |
| \(x[n-k]\) | \(z^{-k} X(z)\) |
| \(y[n-k]\) | \(z^{-k} Y(z)\) |
Given:
\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]
Apply \(\mathcal{Z} \{ \cdot \}\):
\[Y(z) - a_1 z^{-1} Y(z) - a_2 z^{-2} Y(z) = b_0 X(z) + b_1 z^{-1} X(z)\]
Group:
\[Y(z)(1 - a_1 z^{-1} - a_2 z^{-2}) = X(z)(b_0 + b_1 z^{-1})\]
Divide both sides:
\[H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1}}{1 - a_1 z^{-1} - a_2 z^{-2}}\]
Given:
\[y[n] - 0.9 y[n-1] = x[n] - 0.5 x[n-1]\]
Z-Transform:
\[Y(z)(1 - 0.9 z^{-1}) = X(z)(1 - 0.5 z^{-1})\]
Transfer Function:
\[H(z) = \frac{1 - 0.5 z^{-1}}{1 - 0.9 z^{-1}}\]
Let’s analyze \(H(z)\):
Pole-Zero Plot
Visualizes system behavior
Check for: - Stability (poles inside unit circle) - Frequency shaping
Convert this equation:
\[y[n] = 0.6 y[n-1] + x[n] + x[n-1]\]
Find: - \(H(z)\) - Poles and zeros - Plot them in the Z-plane
Notch Filter
A common issue in biosignal processing is removing power‐line interference (50/60 Hz) from, for example, EEG or ECG signals. A simple digital filter to eliminate 60 Hz interference (assuming a sampling frequency \(f_s = 5000\) Hz) is to place complex‐conjugate zeros at
\[ z = e^{\pm j 2\pi\frac{60}{5000}}. \]
The resulting transfer function can be written as
\[ H(z) = 1 \;-\; 2\cos\!\Bigl(2\pi\frac{60}{5000}\Bigr)\,z^{-1} \;+\; z^{-2}. \]
This \(H(z)\) has zeros at \(e^{\pm j2\pi(60/5000)}\) that precisely cancel the 60 Hz component, thereby implementing a notch filter. Moreover, it is a second‐order FIR filter with symmetric coefficients, which grants it linear phase (important to avoid waveform distortion).
def design_filter(zeros=None, poles=None, gain=1.0):
"""
Diseña un filtro digital a partir de ceros y/o polos y una ganancia.
Parámetros:
- zeros: lista de ceros (raíces del numerador), o None para no incluir
- poles: lista de polos (raíces del denominador), o None para no incluir
- gain: ganancia escalar del filtro
Devuelve:
- b: coeficientes del numerador
- a: coeficientes del denominador
"""
# Si no se pasan ceros, asumimos un FIR trivial (b = [gain])
if zeros:
b = gain * np.poly(zeros)
else:
b = np.array([gain], dtype=float)
# Si no se pasan polos, asumimos sistema FIR (a = [1])
if poles:
a = np.poly(poles)
else:
a = np.array([1.0], dtype=float)
return b, adef simulate_ecg(duration=10.0, fs=500, heart_rate=60):
"""
Simula un ECG sintético basado en la superposición de ondas gaussianas.
Parámetros:
- duration: duración de la señal en segundos
- fs: frecuencia de muestreo en Hz
- heart_rate: frecuencia cardiaca en latidos por minuto
Devuelve:
- t: vector de tiempos
- ecg: señal simulada de ECG en mV
"""
dt = 1 / fs
t = np.arange(0, duration, dt)
rr = 60 / heart_rate # intervalo RR en segundos
# Inicializar señal
ecg = np.zeros_like(t)
# Parámetros de las ondas (posiciones relativas en segundos)
# P wave
p_amp, p_dur, p_delay = 0.25, 0.09, 0.16
# Q wave
q_amp, q_dur, q_delay = -0.05, 0.066, 0.166
# R wave
r_amp, r_dur, r_delay = 1.6, 0.1, 0.166
# S wave
s_amp, s_dur, s_delay = -0.25, 0.066, 0.19
# T wave
t_amp, t_dur, t_delay = 0.35, 0.142, 0.36
# Generar cada latido
for beat_start in np.arange(0, duration, rr):
mask = (t >= beat_start) & (t < beat_start + rr)
tb = t[mask] - beat_start
ecg[mask] += gaussian(tb, p_delay, p_dur / 2, p_amp)
ecg[mask] += gaussian(tb, q_delay, q_dur / 2, q_amp)
ecg[mask] += gaussian(tb, r_delay, r_dur / 2, r_amp)
ecg[mask] += gaussian(tb, s_delay, s_dur / 2, s_amp)
ecg[mask] += gaussian(tb, t_delay, t_dur / 2, t_amp)
return t, ecg+0.1*np.cos(2*np.pi*60*t)# Parámetros de simulación
DURATION = 10 # segundos
FS = 500 # Hz
HR = 70 # latidos por minuto
# Generar señal
t, ecg_signal = simulate_ecg(duration=DURATION, fs=FS, heart_rate=HR)
# Graficar resultado
plt.figure(figsize=(12, 4))
plt.plot(t, ecg_signal, linewidth=1)
plt.title(f'Señal de ECG sintética ({HR} bpm)')
plt.xlabel('Tiempo (s)')
plt.ylabel('Amplitud (mV)')
plt.grid(True)
plt.tight_layout()
plt.show()FIR filters
FIR filters (Finite Impulse Response) are widely used in biomedical processing because they can be designed to have linear phase response, avoiding phase distortion in the filtered signal (which is useful for preserving the morphology of ECG waves, for example).
Filter Design Process
Defining the desired ideal frequency response \(H_d(e^{j\omega})\).
Obtaining the ideal impulse response \(h_d[n]\) as the inverse Fourier transform of \(H_d\).
Truncating \(h_d[n]\) (which is usually infinite or very long) with a window function \(w[n]\) to obtain a realizable FIR filter
\[ h[n] = h_d[n]\,w[n]. \]
You obtain \(h_d[n]\) as the inverse discrete–time Fourier transform of your ideal frequency response \(H_d(e^{j\omega})\). For a low-pass filter with cutoff \(\omega_c\),
\[ H_d(e^{j\omega}) = \begin{cases} 1, & |\omega|\le\omega_c,\\ 0, & \text{otherwise}. \end{cases} \]
By definition of the inverse DTFT,
\[ h_d[n] \;=\; \frac{1}{2\pi}\int_{-\pi}^{\pi} H_d(e^{j\omega})\,e^{j\omega n}\,d\omega \;=\;\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} e^{j\omega n}\,d\omega. \]
Carry out the integral in two cases:
For \(n\neq 0\):
\[ h_d[n] =\frac{1}{2\pi}\,\frac{e^{j\omega n}}{j\,n}\Biggr|_{-\omega_c}^{\omega_c} =\frac{1}{2\pi}\,\frac{e^{j\omega_c n}-e^{-j\omega_c n}}{j\,n} =\frac{\sin(\omega_c\,n)}{\pi\,n}. \]
For \(n=0\):
\[ h_d[0] =\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} 1\,d\omega =\frac{2\,\omega_c}{2\pi} =\frac{\omega_c}{\pi}. \]
Putting both together,
\[ h_d[n] =\begin{cases} \dfrac{\sin(\omega_c\,n)}{\pi\,n}, & n\neq 0,\\[1em] \dfrac{\omega_c}{\pi}, & n=0. \end{cases} \]
If your cutoff is specified in Hz, \(f_c\), and sampling rate is \(f_s\), then
\[ \omega_c = 2\pi\,\frac{f_c}{f_s}, \]
so you can write
\[ h_d[n] =\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{\pi\,n} =\;2\;\frac{f_c}{f_s}\;\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{2\pi \frac{f_c}{f_s}\,n} =2\frac{f_c}{f_s}\,\mathrm{sinc}\!\Bigl(2\frac{f_c}{f_s}\,n\Bigr), \]
where we define the normalized sinc as \(\mathrm{sinc}(x)=\frac{\sin(\pi x)}{\pi x}\).
The typical characteristics of classic windows are:
Rectangular: narrowest main lobe (width ≈ 4π/M rad) but highest sidelobes (first sidelobe ≈ –13 dB, stop-band attenuation ≈ 21 dB). It gives the steepest transition for a given filter order, at the expense of poorer stop-band rejection.
Bartlett (triangular): somewhat wider main lobe (≈ 8π/M), sidelobes ≈ –25 dB.
Hann: main lobe ≈ 8π/M, sidelobes ≈ –31 dB (better rejection than rectangular but smoother transitions).
Hamming: main lobe ≈ 8π/M, sidelobes ≈ –41 dB (minimum stop-band attenuation ≈ 53 dB). Very popular for its good compromise between transition width and stop-band attenuation.
Blackman: wider main lobe (≈ 12π/M) but very low sidelobes (≈ –57 dB, attenuation ≈ 74 dB).
Kaiser: allows selection of a parameter β to control sidelobe attenuation, offering flexibility. Approximately, to achieve A dB of attenuation,
\[ \beta \approx 0.1102\,(A - 8.7)\quad(\text{for }A>50), \]
and the normalized transition width Δω relates to the order M and β by
\[ M \approx \frac{A - 8}{2.285\,\Delta\omega} \]
These formulas stem from Kaiser’s approximations and help in sizing the filter.
(-100.0, 5.0)
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import (
boxcar, # Rectangular
bartlett, # Triangular
hann, # Hann
hamming, # Hamming
blackman, # Blackman
kaiser # Kaiser
)
# Parameters
M = 64 # window length
beta = 8.6 # Kaiser parameter for moderate sidelobe attenuation
nfft = 512 # FFT length for frequency response
fs = 1.0 # normalized sampling rate
# Generate windows
windows = {
'Rectangular': boxcar(M),
'Bartlett': bartlett(M),
'Hann': hann(M),
'Hamming': hamming(M),
'Blackman': blackman(M),
f'Kaiser (β={beta})': kaiser(M, beta)
}
# Time-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
plt.plot(np.arange(M), w, label=name)
plt.title('Window Functions — Time Domain')
plt.xlabel('Sample index n')
plt.ylabel('Amplitude')
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()
# Frequency-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
# Compute normalized frequency response
W = np.fft.fft(w, nfft)
W_mag = 20 * np.log10(np.abs(W) / np.max(np.abs(W)))
freqs = np.fft.fftfreq(nfft, d=1/fs)
# Only plot 0 ≤ f ≤ 0.5 (normalized Nyquist)
idx = np.logical_and(freqs >= 0, freqs <= 0.5)
plt.plot(freqs[idx], W_mag[idx], label=name)
plt.title('Window Functions — Frequency Response')
plt.xlabel('Normalized Frequency (cycles/sample)')
plt.ylabel('Magnitude (dB)')
plt.ylim(-100, 5)
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()
plt.show()N = 101
fc = 30
# 2. Compute the ideal impulse response h_d[n] of a low-pass filter
n = np.arange(N)
M = (N - 1) / 2
# Using the normalized sinc: sinc(x) = sin(pi*x)/(pi*x)
h_d = (2 * fc / FS) * np.sinc(2 * fc * (n - M) / FS)
# 3. Choose a window w[n] (here: Hamming) and truncate
w = np.hamming(N)
h = h_d * w
# 4. Normalize to ensure unity gain at DC
h /= np.sum(h)
# 6. Filter it (only allowed library call)
y = sig.lfilter(h, 1.0, ecg_signal)
# 7. Plot input vs. output
plt.figure(figsize=(8, 4))
plt.plot(t, ecg_signal, label='Original signal')
plt.plot(t, y, label='Filtered signal', linewidth=2)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Low-pass FIR (window method) – Cutoff 100 Hz')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()